#include <math.h>
#include <stdio.h>

#ifndef M_PI
#define M_PI   3.14159265358979323846
#endif 

/**
 * @note: 黄道的严格定义是: 地月系质心绕太阳公转的瞬时平均轨道平面与天球相交的大圆. 黄道是天球上黄道坐标系的基圈.
 * @note: 太阳的赤纬等于太阳入射光与地球赤道之间的角度
 *        太阳赤纬的最高度数为23°27', 夏至时太阳的赤纬为+23°27',
 *        冬至时太阳的赤纬为-23°27'. 春分和秋分时太阳的赤纬为0°.
 */

double sin_degree(double degree)
{
    return sin(degree / 360.0 * 2 * M_PI);
}

double cos_degree(double degree)
{
    return cos(degree / 360.0 * 2 * M_PI);
}

double asin_degree(double value)
{
    return asin(value) * 360.0 / (2 * M_PI);
}

double acos_degree(double value)
{
    return acos(value) * 360.0 / (2 * M_PI);
}

double getCurrentJDate(int year, int month, int day,
                       int hour, int minute, int second)
{
    double a = floor((14 - month) / 12);
    double y = year + 4800 - a;
    double m = month + 12 * a - 3;

    return day + floor((153 * m + 2) / 5) + 365 * y +
           floor(y / 4) - floor(y / 100) + floor(y / 400) - 32045 - 0.5 +
           (double)hour / 24 + (double)minute / 24 * 60 +
           (double)second / 24 * 60 * 60;
}

int main(int argc, char *argv[])
{
    //经度(西经为正)
    double lw = - 121.0 + 29.0 / 60.0;
    //纬度(北纬为正)
    double phi = 31.0 + 11.0 / 60.0;
    //当前儒略日
    double jdate = getCurrentJDate(2011, 11, 30, 12, 0, 0);
    //n为以2000年1月1日为基准(0)的儒略世纪数
    double n = round(jdate - 2451545 - 0.0009 - lw / 360);
    //jate_star是地面经度的太阳正午大约时刻
    double jdate_star = 2451545 + 0.0009 + lw / 360 + n;
    //太阳的平近点角
    double m = fmod(357.5291 + 0.98560028 * (jdate_star - 2451545), 360.0);
    //中心差
    double c = 1.9148 * sin_degree(m) + 0.0200 * sin_degree(2 * m) + 0.0003 * sin_degree(3 * m);
    //黄经
    double lambda = fmod(m + 102.9372 + c + 180, 360);
    //过中天
    double jdate_transit = jdate_star + (0.0053 * sin_degree(m)) - (0.0069 * sin_degree(2 * lambda));
    //太阳赤纬
    double delta = asin_degree(sin_degree(lambda) * sin_degree(23.45));
    //时角
    double w_0 = acos_degree((sin_degree(-0.83) - sin_degree(phi) * sin_degree(delta)) / (cos_degree(phi) * cos_degree(delta)));
    //日落
    double jdate_set = 2451545 + 0.0009 + ((w_0 + lw) / 360.0 + n + 0.00053 * sin_degree(m)) -
                       0.0069 * sin_degree(2 * lambda);
    //日出
    double jdate_rise = jdate_transit - (jdate_set - jdate_transit);

    printf("w_0: %f\n", w_0);
    printf("jdate_rise: %f\n", jdate_rise);
    printf("jdate_set: %f\n", jdate_set);

    return 0;
}
